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Ergodic properties of a generic non-integrable quantum many-body system in 

thermodynamic limit 

Tomaz Prosen 

Physics Department, Faculty of Mathematics and Physics, University of Ljubljana, Jadranska 19, 1111 Ljubljana, Slovenia 

(February 1, 2008) 

We study a generic but simple non-integrable quantum many-body system of locally interacting par- 
ticles, namely a kicked t — V model of spinless fermions on 1-dim lattice (equivalent to a kicked 
Heisenberg XX- Z chain of 1/2 spins). Statistical properties of dynamics (quantum ergodicity and 
quantum mixing) and the nature of quantum transport in thermodynamic limit are considered as 
the kick parameters (which control the degree of non-integrability) are varied. We find and demon- 
strate ballistic transport and non-ergodic, non-mixing dynamics (implying infinite conductivity at 
all temperatures) in the integrable regime of zero or very small kick parameters, and more generally 
and important, also in non-integrable regime of intermediate values of kicked parameters, whereas 
only for sufficiently large kick parameters we recover quantum ergodicity and mixing implying nor- 
mal (diffusive) transport. We propose an order parameter (charge stiffness D) which controls the 
phase transition from non-mixing/non-ergodic dynamics (ordered phase, D > 0) to mixing/ergodic 
dynamics (disordered phase, D — 0) in the thermodynamic limit. Furthermore, we find exponential 
decay of time-correlation function in the regime of mixing dynamics. 

The results are obtained consistently within three different numerical and analytical approaches: 
(i) time evolution of a finite system and direct computation of time correlation functions, (ii) full 
diagonalization of finite systems and statistical analysis of stationary data, and (iii) algebraic con- 
struction of quantum invariants of motion of an infinite system, in particular the time averaged 
observables. 

PACS numbers: 05.45. +b, 05.30.Fk, 72.10.Bg 
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I. INTRODUCTION 

It has been a common belief for a long time, that a 
large system of sufficiently many interacting particles 
should fill uniformly the entire available phase space. 
This is known as ergodic hypothesis, one of the corner- 
stones of the statistical mechanics, and is a necessary 
assumption to justify the use of canonical ensembles and 
derivation of fundamental laws of statistical physics, such 
as transport laws (e.g. Ohm's law or Fourier's law). 

However, the proof together with the precise condi- 
tions for the validity of ergodic hypothesis is still one 
of the most fundamental unsolved problems of theoreti- 
cal physics. Even in the context of purely classical dy- 
namics, the ergodic theory pL0|, though it is an involved 
and beautiful mathematical discipline, can make strong 
statements only for a very limited class of systems, while 
generic dynamical systems, especially those consisting of 
many interacting particles, are far from being understood 
Even less is known about ergodic properties of 
generic quantum many body systems, which is precisely 
the objective of this paper. A closed (finite) and bounded 
quantum system of size L and with a finite number N of 
particles has a discrete spectrum, hence its time evolution 
is quasi-periodic, and accordingly it is non-ergodic and 
non-mixing, as we shall define below. However, in the 
thermodynamic limit (TL), of diverging size L — > oo and 
density of particles p = N / L fixed, the spectrum of the 
quantum propagator may accuire a continuous compo- 



nent, and one may expect genuine properties of quantum 
ergodicity and quantum mixing to set in provided the 
strength of non-linear interaction is sufficiently strong. In 
this paper we deal with general non-autonomous many- 
body systems with Hamiltonians H{t) which explicitly 
depend on time r. Therefore the entire Hilbert space of 
many-body quantum configurations (Fock space) is dy- 
namically accessible and the 'micro-canonical' average of 
an intensive observable represented by an operator A, 
reads 



(A) = lim 

L >cx 



trA 
tr 1 ' 



(1) 



If the system possesses a group of exact geometric or 
dynamical symmetries, the trace in eq. (Q) may be con- 
sidered only over a specific symmetry class of the Fock 
space w.r.t. symmetry group. For example, if the system 
is autonomous, energy is conserved, and eq. (Q) should 
be replaced by the average over a specific "energy shell" , 
or, as often, if the number N of particles (or the particle 
density p — N/L) is preserved, then the micro-canonical 
average should be performed over the Fock supscace of 
N-particle configurations 



lim 



tr(A5 [pL]:A r) 



trS 



[pL],N 



(2) 



where [x] is an integer part of x and 8 m ,n is a Kronecker 
symbol. When we will like to keep the size L in the aver- 
age @ fixed and finite we would write (A) p. Although 
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in this abstract discussion we would like to avoid the no- 
tion of temperature ||, one may think of ([l]) or (|^) as 
canonical averages at very large or infinite temperature, 
(3 = {ksT)- 1 - 0. 

The system is quantum ergodic if the time average of 
an arbitrary observable in Heisenberg picture A(t) equals 
to the micro-canonical average (A) times a unit operator 
(over the corresponding desymmetrized Fock subspace) 



A := lim — / drAfr) 



(A)l. 



(3) 



In case where one has a constant of motion, e.g. density p 
(or energy E = H, etc.) one should define the ergodicity 
through the spectral resolution of the relevant invariant 
operator, p = J p'dE p > , namely 



A = J (A) p ,dE p 



(4) 



When the micro-canonical average does not depend on 
the eigenvalues of the symmetry operations or quantum 
numbers, e.g. (A) p — (A), definition (Q) is equivalent to 
a simple one (g). 

Even stronger ergodic property is quantum mixing, 
which is defined very generally || as follows: The infinite 
(L = oo) quantum many-body system is called quantum 
mixing if time-correlations for an arbitrary pair of quan- 
tum observables in Heisenberg representation, A(t) and 
B(t), decay to zero 



C AB (r) := (A(r)B(O)) 
lim Cab (t) = 0. 



(A)(B), 



(5) 



Note that formula (||) implies that TL should be consid- 
ered prior to the time limit, r — > oo, since these two limits 
do not generally commute ||. Again, in case of addi- 
tional symmetry, mixing over separate symmetry classes 
can be studied; as well as uniform mixing over entire Fock 
space (which makes sense when canonical averages of A 
and B do not depend on quantum numbers like p). 

In Ref. H quantum mixing of a system of interacting 
bosons has been related to a hard chaos of the corre- 
sponding classical (mean field) model. However, general 
quantum systems need not possess the classical limit; 
when they do, the general definitions (J3|j|) go over to 
the correct definitions of classical ergodicity and mixing 
of the corresponding classical counterparts. 

It is easy to see that the same implication holds as 
in classical mechanics 0] : Quantum mixing (||) implies 
quantum ergodicity (0)7 To see this, observe that as 
a simple consequence of (||), time averaged correlation 
function should vanish, 



Cab ■= lim — 

T-*oo T 

This fact is equivalent to 



dTC AB {T) =0. 



Then we immediately see that observable in brackets 
should vanish A — (A)l = 0, since observable B is ar- 
bitrary, q.e.d. The last argument can be reversed, so we 
see that quantum ergodicity is equivalent to Cab — 
for an arbitrary pair A, B. 

We expect that quantum mixing implies universal sta- 
tistical properties of energy spectra (and also universal 
statistics of occupation numbers ffl], matrix elements, 
etc) described by Random Matrix Theory |[, (RMT). 
Random matrix spectral statistics have indeed been 
demonstrated numerically for few strongly non-integrable 
many-body systems Q. On the other hand, completely 
integrable quantum many-body systems (having an in- 
finite set of independent conservation laws Q n ,n = 
1,2,3...) are obviously non-ergodic (Q n — Q n ^ (Q) n l), 
and therefore non-mixing, and characterized by the uni- 
versal Poissonian spectral statistics 0L 

It has been pointed out recently |l(],[ll| that integra- 
bility typically implies non- vanishing stiffness, i.e. ideal 
conductance with infinite transport coefficients (or ideal 
insulating state). Indeed, there is a direct implication 
of quantum mixing and quantum ergodicity on quantum 
transport. One should simply inspect a Kubo formula 
]l2]| , which relates the real part of the transport coeffi- 
cient, e.g. electric conductivity & '(u>), to the cosine trans- 
form of the autocorrelation function of the electric cur- 
rent observable J, written for high temperatures (small 
(3) as 



a'(uj) = 1/3 cos(wr)(ij(0)J(r))dT. 



(6) 



The transport is diffusive and the system behaves as a 
normal conductor if zero-frequency (d.c.) conductivity 
is finite, cr'(0) < oo, which means that the time inte- 
gral of the current-current correlation function should be 
finite; this is true if the system is mixing and if time cor- 
relations decay sufficiently fast, e.g. it is sufficient that 
|(±J(0)J(t))| < C\r\- a for some C > and a > 1. 

On the other hand, if the transport is ballistic, d.c. 
conductivity diverges cr'(0) = oo and frequency depen- 
dent conductivity can be written as a sum of delta-spike 
and a regularized conductivity 



er'(cj) = D5(ui) + <y reg (w), 



D = lim 

e-+0 



a'{u)du (7) 



The weight D is known as a charge stiffness (or Drude 
weight) and is proportional to the averaged current- 
current time correlator 



D = 0Dj, D A 



lim lim 

T^oo L^oo 2TL 



C L AA (j)dr (8) 



— T 



(AB)-(A)(B) = ((A-(A)1)B) = 0, 



Therefore, non- vanishing charge stiffness Dj ^ 0, mean- 
ing a ballistic electric transport, is a sufficient condition 
for deviation from quantum ergodicity, so Dj will be 
extensively used as a quantitative indicator of quantum 
(non)crgodicity throughout the rest of this paper. 
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In a generic integrable system one can find an (infi- 
nite) set of invariant extensive observables, the so-called 
conserved charges Q n . Mazur and Suzuki ||3 have 
proposed a 'Parseval-like' inequality for the time aver- 
aged autocorrelator of any extensive observable A 



(yAQ n ) 



(9) 



using any suitable (sub)set of conserved charges {Q m }, 
such that {j^QnQm) — if n 7^ m. For an integrable 
system one therefore proves ideal (ballistic) transport 
(Dj > 0) if at least one term in an infinite sum on RHS 
of eq. (||) applied to the current observable A = J is 
non- vanishing (what is typically the case) . It is con- 
venient to say that {Q n } is a complete set of conserved 
charges if (^) is an exact equality for any observable A 
(which is 'square-summable', (j^A 2 ) < 00). 

The important and delicate question is whether non- 
crgodicity of an integrable system in TL can be struc- 
turally stable against generic and finite non-integrable 
perturbation. In this paper we will present clear numeri- 
cal evidences based on various different and independent 
numerical methods in support of a conjecture claiming 
an affirmative answer to the above question. 

Conjecture: Let H\ be a continuous family of generic 
infinite quantum-many body systems with local interac- 
tion , such that Hq is completely-integrable while H\ 
are non-integrable for almost any A 7^ 0. Then, 3A C , such 
that H\ are non-ergodic and non-mixing for |A| < A c . 

One important consequence of the above Conjecture 
would be a large class of not only completely integrable 
but also nearly-integrable many-body systems for which 
ideal transport and infinite conductance would be ex- 
pected. 

In Sec. II we define a two parametric family of generic 
non-integrable many-body systems pi , namely a kicked 
t-V model of interacting spinless fermions, or equiva- 
lently, a kicked Heisenberg XX-Z spin-i chain, and de- 
scribe the basic properties of the model. In Sec. Ill we 
show how efficient explicit time evolution of the above 
model of finite but quite large size L can be computed, 
and present results on extensive numerical computation 
of time correlation functions. By letting the size L 
to increase and inspecting TL, we clearly identify two 
regimes of quantum motion: non-mixing regime for small 
and intermediate values of kick parameters where time- 
correlation functions typically saturate to constant non- 
vanishing values, and exponentially mixing regime for suf- 
ficiently large values of kick parameters where time cor- 
relation functions decay exponentially. As a complemen- 
tary approach to direct time evolution (time domain) we 
perform, in Sec. IV, complete diagonalization of station- 
ary problem (frequency domain) for finite sizes L. Using 
stationary data we have computed and analyzed short- 
range and long-range quasi-energy level statistics. By 



means of matrix elements of the current observable J 
we have also directly computed the conductance and the 
charge stiffness Dj. The obtained results are in quanti- 
tative agreement with a direct time evolution (Sec. III). 
In Sec. V we outline completely different and inde- 
pendent method of computing time-averaged observables 
and quantitative indicators of quantum ergodicity such as 
the charge stiffness Dj, by making use of extensive com- 
puterized Lie algebra. This third method, in contrast to 
other two, refers directly to infinite systems (infinite lat- 
tices L = 00) and, again, gives compatible results. In 
Sec. VI we give more arguments in support of our Con- 
jecture by discussing relevant published and unpublished 
results and conclude. 



II. THE MODEL 

No general analytical methods exist to deal with dy- 
namics of non-integrable quantum many-body systems. 
A beautiful solition theory based on inverse scattering 
and algebraic Bethe ansatz |l7| is unfortunately applica- 
ble only to a limited class of very special, namely com- 
pletely integrable many-body systems. Therefore one is 
left to numerical experiment to learn about dynamics of 
generic quantum systems, encouraged with the fact that 
numerical and experimental investigations of dynamical 
systems of one or few particles has been a very fruitful 
area of research (known as Quantum Chaos) over the past 
twenty years pl|. However, one should be very careful 
in picking out the toy model, since the fact, that the di- 
mensionally of the Hilbert space (Fock space of quantum 
many-body states) grows exponentially with increasing 
system size L, makes a serious quantitative study of TL 
almost prohibitive. Here we propose the simplest many- 
body system that we can think of: one dimensional lat- 
tice of spinless fermions of size L, and for the reasons 
which will become clear in the next section, we decide 
to break integrability by taking time-dependent interac- 
tion which is switched on periodically by means of 5- 
kicks. Therefore, the time-dependent Hamiltonian of our 
"kicked t-V model" (KtV) fy§ reads 

L-l 

H(t) = [-2*( c i c J-+l + h - c -) + W Vn j n j+l] ■ ( 10 ) 

3=0 

Cj, Cj are fermionic creation and annihilation op- 



erators satisfying canonical anticommutation relations 

[ c j, c k} + ■= CjCk + C k Cj = 0, [ct,Cfe] + = 6jk, Tlj = c]cj 

are number operators, and periodic boundary condi- 
tions are imposed c L = c . S p (t) = J2m=-oc^( T ~ m ) 
is a periodic <5-function. We use units in which h — 
(time between collisions) = (lattice spacing) = 1. The 
hopping amplitude t and interaction strength V are in- 
dependent (kick) parameters. An important and useful 
property of kicked systems like (|l0|) is the fact that the 
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evolution propagator over one period (Floquet operator) 
factorizes into the product of kinetic and potential part 

= exp(-iW) exp(-iT), (11) 

L-l 

'c j+ i), (12) 



3=0 
L-l 

3=0 



(13) 



We have used a Pierels phase <j) in order to introduce a 
particle current 



L-l 



j = (i/t)u\d/d4>)u\\<t>=o = (4 c i+i - 4+1^) 

(14) 



fc=0 



(which is divided by t for convenience); elsewhere we put 
<fi := 0. Note that the kinetic energy T as well as the 
current J are diagonal in momentum representation 



L-l 

T = t^(l- cos(sfc))n fe , 

fc=0 
L-l 

J = sin(sfc)nfe, 

Ai=0 



(15) 
(16) 



where s = 2tt/L, while tilde refers to momentum repre- 
sentation of field operators 



L-l 



c k = L 1/2 ^ exp(isjk)cj, h k = c\c k . (17) 



3=0 



Using a well know Jordan- Wigner transformation ]l9| ] 
one can map 1-dim lattice of spinless fermions to a 
spin— | chain described by Pauli operators — (crj ± 
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f)/V% and a?, namely 



j'-i 



= V%Cj exp I in rij' 

2'=0 



fj| = 2nj - 1, 



(18) 



which satisfy canonical commutation relations 

Wj> a l\ = i £ e»vr,e]Sjk, V>>v,r)£ {x, y, z}. 



In fact, Jordan- Wigner transformation (|I|) maps KtV 
model on a kicked Heisenberg XX-Z chain 



H(t)=T + 6 p (t)W, 



L-l 



J=0 



L-l 



^ = ^>>X + i + K)- 



4 V"j "jH 



(19) 
(20) 

(21) 



The last term of potential ( pl| ) is irrelevant since the total 
z— spin S z — J2f=o a j — ^ — \L is a constant of motion, 
[U,S z ]=0. 

Interaction strength V is cyclic parameter V = V 
(mod 2tt), since the spectrum of W/V is a set of inte- 
gers (see eq. (p"3|)). KtV model is integrable/solvable in 
three special (limiting) cases: 

1. t = Q, 1-dim Ising model, 

2. V = (mod 27r), 1-dim free fermions, or equiva- 
lently, 1-dim Heisenberg XX ^-spin chain, 

3. tV — > and A = t/V finite, (continuous-time) 1- 
dim Heisenberg XXZ i-spin chain. 

For t ^ 0, V 7^ (mod 2tt), KtV is expected to be non- 
integrable, possibly quantum ergodic and mixing. 

To conclude this section, let us list the symmetries of 
a general KtV model (for arbitrary t and V): In addition 
to the trivial conservation law, namely the number or 
density of particles 



N 



L-l 

£ n J. 

3=0 



N 

p=-, [U,N]=0, (22) 



and the total quasi-momentum K G {0,1,..., L — 1} 
which is defined as an eigenvalue of a unitary transla- 
tional symmetry operation S 



I L-l N 

S = cxp(isK) = exp I is kn k 

\ k=0 / 



[U,S]=0, (23) 



the KtV model has two (geometric) 'reflection' symme- 
tries, namely the parity 

V:r, •(•/ ,. VU = UV, V 2 = l, (24) 

and, for even size L, the particle-hole transformation 



1Z : Cj 



(-iy c 



KU = UK, 1Z 



1. 



(25) 



Note on notation: Symbols wearing a hat " denote lin- 
ear transformations over the operator space of quantum 
observables. 



III. FIRST METHOD: DIRECT TIME 
EVOLUTION AND CORRELATION FUNCTIONS 

For a fixed size L and fixed number of fermions N, 
a unitary quantum many-body map U ( p"l| ) acts over a 
Fock space of dimension 
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LI 



L 



Nj Nl(L-N)\ 



(26) 



Dynamics of a given initial many-body state IV'(O)) is 
given by a simple iteration of the Floquet map 



\if){m)) = U\^{m - 1)) = U m \ip(0)). 



(27) 



Many body states can be expanded in a complete 
basis of the Fock space (of Slater determinants), for which 
we may choose either the position states, labelled by sets 
of N ordered integers j = (ji,... , jjv) 

|J> = % ■ • • c jN \0), < j x < . . . < j N < L, (28) 

or the momentum states, labeled by sets of N ordered 
integers k 

\k) =c kl ---h N \0), < fci < ... < k N < L. (29) 

An important observation, implicitly made already in 
previous section (|T^Jl^,|l^) , is that kinetic propagator 
cxp(— iT) is diagonal in momentum representation while 
potential propagator exp(— iW) is diagonal in position 
representation 



If-, := (fc|exp(-zT)|fc') 



N 



H,k' ex P it^2,{cos(sk n ) - 1) 



n=l 



U&:= (j\eM-iW)\f) 



N 



n,n'=l 



(31) 



Therefore, one may formulate a very efficient algorithm 
to perform explicit time evolution of many-body states, 
provided that it is possible to switch between the two 
representations ( |28| ) and ( p9| ) as efficiently as in the prob- 
lem of one kicked quantum particle N — 1, e.g. kicked 
rotor |20) or kicked Harper model j2l[] , namely by means 
of Fast Fourier Transformation (FFT) algorithm. In- 
deed, we succeeded to develop a fast algorithm which 
performs such anti-symmetrized multi-dimensional dis- 
crete Fourier transformation 



6\k) 



(32) 



in roughly J\f\og 2 J\f Floating Point Operations (FPO). 
It is based on a factorization of L— site discrete Fourier 
transformation into the product of ~ L log 2 L 2-site 
transformations parametrized with 2x2 sub-matrices 
(a,/3;j,S)jf, which are successively applied to creation 
operators, (cj,ct) <— (act + (3cj, , jCj + 6c^, ) , in all Slater 

determinants IT ra cj n |0) which contain a particle at sites 
j or j' . (One should be careful in dealing with fermionic 
signs of Slater determinants when sorting the factors in 



the product n„cjjO)). In case when L — 2 P , factoriza- 
tion of FFT to a chain of 2-site (in such case unitary) 
transformations is easily deduced by inspecting a con- 
ventional FFT algorithm (such as the one implemented 
in |^2|), while for more general lattice sizes (we have so 
far implemented our Fermionic FFT (FFFT) algorithm 
also for L = 10, 12, 15, 20, 24, 30, 40) we factorized the op- 
timal schemes developed by Winograd Q. Our FFFT 
algorithm requires almost no extra storage apart from a 
vector of J\f c- numbers where quantum many-body state 
is stored. Therefore, the map ( ll]) is iterated on a vector 
-0g(m) = (fc|-0(m)), using the matrix composition 



U 



F*U W FU 7 



(33) 



in roughly 2A/" log 2 TV" FPO per time step which is by 
far superior to 'brute-force' methods based on complete 
diagonalization of U and expansion of time-evolving state 
\ip(m)} in terms of eigenstates of U. 

Let us now consider time autocorrelation functions 
of two 'generic' observables, namely the current J and 
rescaled traceless kinetic energy T 1 



C L j(m) 



(J(0)J(m)) L p 



(34) 



CHm) := -(T'(0)T'(m)) L p , T' = -T - N 



L 



(30) J( m ) = W m JU m , T'(m) = W m T'U m . Note that 
(J)p = (T')p = 0. These two observables belong to 
different symmetry classes with respect to the parity op- 
eration 



VJ = -JV, VT' = T'V, 



(35) 



so we choose both to check whether many-body dynam- 
ics can substantially depend of the symmetry class w.r.t. 
'reflection' symmetry. Conveniently, both observables, J 
and T' , are diagonal in momentum representation, 



J\k) = J % \k), Jg=^sin(sfc n ) I (36) 

ra=l 
JV 

T'\k)=T' % \k), Tt = -J2™s(sk n ), (37) 



so the time auto-correlation functions can be computed 
from time evolution of a (complete) set of TV 7 (= TV) 
initial momentum states k! 

c ' A{m) = Ljpil A ^T. A iPrA m ) (38) 



where A is any observable which is diagonal in momen- 
tum basis (here either J or T") and 



p m ,{rn) = \{mm))\ 2 = \{k\U m \k')f 



(39) 
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When dimensionality Af becomes prohibitively large we 
suggest to estimate the microcanonical averages (p4) by 
taking a smaller 1 <C Af' <C Af but uniformly random 
sample of Af' initial states \k'}, (.) = -fr Epv^l-I^'H 1 + 
0(l/\f Af')). Therefore, numerical computation of corre- 
lation function C'a (m) for m = 1, . . . , M can be per- 
formed in - (2MAf'Af\og 2 Af)/L FPO. Reduction for a 
factor 1/L w.r.t. naive FPO count is due to translational 
symmetry ( |23| ) , since one can simultaneously simulate the 
dynamics of L different states with different values of the 
conserved total momentum K = Yln=i ( m od L). 

Let us for the time being fix the density of particles 
p = N/L = 1/4. We have performed extensive nu- 
merical computations of time-correlation functions (by 
means of explicit time evolution (|3§1)), for sizes L = 
8, 12, 16, 20, 24, 32 (at L = 32 the dimensionality of Fock 
subspace is Af — 10518300), and systematically scanned 
the parameter space (t, V). We have clearly identified 
two regimes where we were able to probe the TL, i.e. 
where time-correlation functions turned out to be stable 
against the variation of the system size L: 

(i) Quantum ergodic and mixing regime for sufficiently 
large value of parameter t and for any value of param- 
eter V (away from 'integrable axis' V — (mod 27r)). 
In this regime, time correlation functions are rapidly de- 
creasing (Case t = V = 4 is shown in Fig.l). However, 
for a finite size L the quantum system is almost never 
mixing, so correlation functions saturate, on a time scale 
p(L), to a small but finite value of the stiffness 



U A 



lim 

M^oc 



1 



2M + 1 



M 
i=-M 



(40) 



(Here A = J or T' .) In order to avoid transient be- 
havior at small times m and incorporating the time re- 
versal symmetry, C\(m) = C^(— m), the time aver- 
ages like ( ^p| ) have been numerically estimated as D\ — 
W+tT^Lm' C\{m). The sufficiently large averaging 
time scale {M' . . . 2M '} = {30 . . . 60} for all sizes L < 32 
(and for most values of parameters t, V) has been deter- 
mined by direct inspection of the correlation functions 
(see Fig.l). In Fig. 2 we plot correlation functions Cj(m) 
for t = V = 4 in semi-log scale, and show that, as the 
size L is increased, the saturation (Thouless) time scale 
p(L) increases, roughly as p(L) ~ L. So, the Thou- 
less time p(L) clearly diverges in TL. Furthermore, Fig. 2 
gives clear numerical evidence (further supported by re- 
sults shown in Figs. 4, 5) of the exponential decay of time 
correlation functions in TL (or for times smaller than 
p(L) in a system of finite size L) 



Ca(iti) cx exp(— A^m), m 1. 



(41) 



Henceforth, the stiffness should also vanish exponentially 
D\ ~ exp(— XaL) as one approaches TL (Such behav- 
ior has been observed also in Ref. indicating ex- 



ponential mixing in the system sudied there |10| ). In- 
deed, in Fig. 3 we examine 1/L scaling of the charge 
stiffness Dj which is (here shown for t = V = 4 and 
t — V = 2) a clear indication of ergodic and mix- 
ing behavior in TL, Df = 0. In this regime, in TL, 
Kubo conductivity a'(uj = 0) = \(5 Ylm=-<x> Cj{m) is 
finite, cr'(0) < 00, and the transport is dissipative. Fur- 
ther, as shown in ]l6| ], time averaged current of ar- 
bitrary initial momentum state \k') averages to zero 
Jp - lim M ^oo(l/M)E™ =1 (fcVM|fc') = 0, and arib- 
trary initial momentum state \k') explores entire acces- 
sible Fock space, i.e. (k\U m \k') are uniformly Gaussian 
pseudo-random numbers when the discrete time m is suf- 
ficiently large, say larger than the quantum mixing time. 
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FIG. 1. Current autocorrelation function Cj{m) against 
discrete time m for quantum ergodic (t = V = 4, lower 
set of curves for various sizes L) and intermediate regime 
(t = V = 1, upper set of curves) with density p =\. Averag- 
ing over entire Fock space is performed, M' = TV, for L < 20, 
whereas random samples of J\f' = 12000, and TV"' = 800 initial 
states have been used for L — 24, and L — 32, respectively. 



(ii) Non-ergodic and Non-mixing regime for parameter 
t ~ 1 (or smaller) and any value of parameter V. Here 
time correlation functions CU(to) do not decay to zero 
but saturate, around a constant typically non-vanishing 
and positive value of the stiffness D\ > (f4(i|), on a 
short time scale which does not depend on the size L (for 
sufficcntly large sizes L). In Fig.l we plot time corre- 
lation functions Cj(m) for t = V = 1. Please observe 
indeed very weak dependence on the size L. In Fig. 3 we 
also show 1/L scaling of charge stiffness Dj for the cases 
t = V = 1 and t = 1 , V = 2 which clearly indicate a 
finite extrapolated (to 1/L = 0) thermodynamic value 
of the stiffness. This should be considered as evidence of 
non-mixing and non-ergodic behavior in TL. Since in this 
parameter ranges KtV model is also non-integrable, we 
sometimes refer to this regime as an intermediate quan- 
tum dynamics. This behavior corresponds to ideal, bal- 
listic transport with infinite Kubo conductivity a = 00. 
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Furthermore, in |16| it has been shown that in this in- 
termediate regime the time averaged (persistent) current 
is non-vanishing and proportional to the initial current 
Je,, Jp — a Jp, ol — 2Dj/[p(l — p)] (which is the most 
direct probe of ideal, ballistic transport), and that an 
arbitrary time-evolving initial momentum state U m \k') 
remain strongly localized in a non-trivial subregion of 
dynamically accessible Fock space. 




5 10 15 20 25 30 35 40 45 50 
m, discr.time 

FIG. 2. The lower set of curves (t = V = 4) of Fig.l 
in semi-log scale. To emphasize the exponential correlation 
decay we plot also the best exponential fit to the tail of 
Cj" 32 (m) (dash-dotted line). 
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FIG. 3. Stiffness Dj vs. 1/L at constant density p = \ 
and for different values of control parameters in quantum mix- 
ing/ergodic, t = V = 4 and t = V = 2, and intermediate, 
i = 1, V = 2 and t — V = I , regime. Other parameters are 
the same as in Fig.f 

One may use a charge stiffness of an infinite system 



Df as an order parameter controlling the dynamical 
phase transition from a disordered phase (quantum er- 
godic/mixing dynamics) characterized by Dy — to an 
ordered phase (non-ergodic/non-mixing dynamics) char- 
acterized by Dy > 0. The transition point is charac- 
terized by diverging correlation time (or mixing time) 
scale, AT 1 , which diverges when one approaches the tran- 
sition from above, say, with parameter t decreasing to- 
wards certain critical curve t c (V). Of course, in the or- 
dered phase, t < t c (V), Dy > 0, the time-correlations 
have infinite range, CU(oo) =/= 0. The transition is il- 
lustrated in Figs. 4, 5 by plotting correlation functions 
for both observables, Cj(m) (Fig. 4) and Cy(m) (Fig. 5), 
for different values of parameter t and fixed parameter 
V = 2. En estimate of the critical parameter here is 
1.4 < t c {V — 2) < 1.5. Observe the exponential decay 
of correlations in all cases, except possibly at very small 
times where other, smaller time scales may become im- 
portant. 
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FIG. 4. Current auto-correlation functions Cj~ 32 (m) for 
different values of parameter t (see legend) and fixed value of 
parameter V . p = 1/4. 
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FIG. 5. Kinetic auto-correlation functions C^ =32 {m) for 
different values of parameter t (see legend) and fixed value of 
parameter V. p = 1/4. Note that the same scale is used than 
in previous Fig. 4. 
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IV. SECOND METHOD: EXACT 
DIAGONALIZATION OF STATIONARY 
PROBLEM OF FINITE SIZE 

In a more complete but brute-force approach one may 
try to exactly diagonalize the matrix of a one-period evo- 
lution (Floquet) propagator U for a finite size L and 
compute interesting dynamical quantities, such as con- 
ductivity a'(ui) and stiffness D\ directly from the spec- 
trum {(fin} and the set of eigenstates {\n}} of KtV map 
U. 



(42) 



U\n) = exp(— itp n )\n), n = 1 . . .TV. 



Again, it is easiest to work in the momentum basis (|2_ 
and to use the translational symmetry to decompose the 
matrix p into blocks (with fixed value of the total 
quasi- momentum K) of dimension Mk ~ M jL. Only for 
blocks with K = and K = L/2 (if L is even) the parity 
24|) commutes with the translation S (23), 



operation V 

and may be used to further reduce the dimensionality of 
irreducible block by factor 2. The matrix Ut r, (for fixed 
K) has been computed, by means of a decomposition j3^ ) 
and FFFT algorithm, in roughly {M/L) 2 \og 2 M FPO, 
and further diagonalized by means of standard routines 
in roughly (7V/L) 3 FPO giving a set of quasi-energies 
{(£>„} and eigenstates (k\n). 



A. Spectral statistics 

In the so called Quantum Chaology of simple (few) 
body non-integrable system there is a famous Conjectue 
of Bohigas, Giannoni and Schmit pi] , supported by nu- 
merous numerical |2E) | and theoretical arguments [p6| , 
claiming that hard chaos (ergodicity, mixing and posi- 
tive Lyapunov exponents) of a classical counterpart re- 
sults in a universal statistical properties of system's 
(quasi)energy spectrum given by the appropriate ensem- 
ble of random matrices B. On the other hand, in- 
tegrable classical classical dynamics results in univer- 
sal Poissonian statistics of (locally) uncorrelated (quasi)- 
energy levels p7[. Intermediate statistics, which are nei- 
ther RMT nor Poissonian, are found [^8|,^9l for systems 
whose classical dynamics is intermediatc(mixed) with 
regular and chaotic motion coexisting in phase space. 
The connection between integrability/non-integrability 
and statistics has recently been investigated also in few 
well known examples of non-linear many-body systems 
(correlated fermions or interacting spin chains) which 
do not possess a well defined classical limit ||. It has 
been shown that quantum integrability, or (strong) non- 
integrability, of the quantum many-body model again 
correspond to Poissonian, or RMT, behavior of level 
statistics, respectively. No attempt has been made there 
however, to understand the intermediate situation, 
or the thermodynamic limit. 



Inspired by Quantum Chaos we have analyzed the sta- 
tistical properties of the quasi-energy spectrum {ip n } of 
KtV model and searched for signatures of ergodicity and 
mixing of the underlying quantum many-body dynamics 
in TL. For comparison with other results, the density will 
be again fixed to p = 1/4 in the numerical presentation 
which follows. 

First, we have analyzed the common short-range statis- 
tic, namely the integrated (cumulative) nearest neighbor 
level spacing distribution, W(S), giving the probability 
that a random normalized spacing between two adjacent 
eigenphases S n — ^j^ifn+i — fn) is smaller than S 



W{S) 



(43) 




FIG. 6. Cumulative quasi-level spacing distributions 
W(S) for quarter-filled p = 1/4 KtV model. In (a) we plot 
W(S) for several different values of parameter t and fixed pa- 
rameter V (see legend) covering the mixing/non-mixing tran- 
sition, and for maximal computable size L = 24. With dot- 
ted curves we plot for comparison the theoretical, Poissonian, 
Semi- Poissonian, and Wigner (COE), distributions. In (b) 
we show W(S) for fixed kick parameters t = 1, V = 2, and 
for different sizes L = 16, 20, 24. In the insets we plot the 
same objects in log-log scale to emphasize the small spacing 
behavior. Please observe the trend towards linear repulsion 
(quadratic for W(S — » 0) oc S 2 ), even in intermediate regime. 



8 



For moderate values of the size L < 20, the average in 
([l3|) has been computed over all M' = N states from all L 
blocks (symmetry classes) of constant quasi-momentum 
K. (Note that blocks for K = K' and K = L - K ' are 
related by a parity transformation ( pi]) and give identi- 
cal (sub)spectra, so one has to diagonalize only [L/2 + 1] 
different blocks of matrix Urp.) However, for larger size 
L = 24 we have already Mr ~ 5608, so we averaged only 
over one class of fixed quasi-momentum, namely K = 1 , 
Af'=M~ Af/L. (We have carefully checked that the 
statistical properties of partial subspectra are indepen- 
dent of the symmetry class labeled by quasi-momentum 
K.) 

In Fig. 6a we show W(S) for size L = 24 and several 
different values of parameter t (and fixed V = 2) cov- 
ering the transition from non-ergodic to ergodic/mixing 
quantum dynamics. We find almost Poissonian behavior 
W P (S) = 1 - exp(-S) for small t and excellent RMT 
behavior Wcoe(S') = 1 — exp(— nS 2 /i) (Wigner surmise 
approximating the statistics of the infinitely-dimensional 
Circular Orthogonal Ensemble COE |8j, due to time- 
reversal symmetry) for t > t c (V). In the (more interest- 
ing) region of intermediate dynamics 1 ~ t < t c (V) we 
find intermediate statistics interpolating between Poisso- 
nian and COE (see Fig. 6). Interestingly, level statistics 
close to the critical point (for t — 1.4, V — 2) seems 
to be well captured by the so-called Semi-Poisson model 
W S p(S) = 1 - (1 + 25) exp(-2S) @ which has been 
recently used to model the critical level statistics of 3D 
Anderson model |u| . Since it is impossible to make state- 
ments about TL of level statistics based on results for a 
fixed size L = 24, we show in Fig. 6b the dependence of 
W(S) on the size L for fixed parameters t = 1,V = 2 
(in the regime of non-ergodic dynamics). Although the 
intermediate W{S) statistic is closer to Poissonian than 
to COE, it moves towards COE as we approach TL (in- 
crease L). This fact eliminates possible fears of accidental 
integrability of KtV model in the claimed intermediate 
regime. 

Second, we have analyzed the long-range spectral 
statistics, namely the number variance 

£ 2 (S) = (n(S) 2 ) -S 2 , (44) 

giving the variance of the number n(S) of normalized 
(unfolded) levels j^^n m a randomly chosen interval 
of length S (Note 'that (n(S)) = S.) For Poissonian 
and COE model we expect Sp(S') = S and Ep OE (»5') w 
(2/-7T 2 ) log^-TrS), respectively. Here, one should note huge 
degeneracies in the integrable limit t = of Ising model 
(which are quite-common in integrable quantum many- 
body models in general). For small kick parameters t, 
we hence find stronger-than-Poissonian level clustering 
causing faster-than-linear growth of S 2 (S) (Fig. 7a). For 
finite size L = 24, strong level clustering affects also long- 
range statistics in the regime with mixing dynamics in 
TL, namely for t = 2.5 we find good agreement with 
COE statistics only for relatively small spectral ranges 



S < >S m ax ~ 10 1 . It has been checked, however, that 
the agreement with COE improves to hold on longer 
quasi-energy ranges (S max increases), as either the kick- 
parameter t or the size L are increased. In the interme- 
diate regime 1 ~ t < t c (V), number variance approaches 
that of an uncorrelated sequence, T, 2 (S) ~ S, as we ap- 
proach TL (see Fig. 7b for case t = 1, V = 2). However, 
for finite L, the phenomenon of saturation sets in p2[ , 
namely when the scaled energy range S — S* is of the 
order of the density of states 

S* = 0.5AA/L, (45) 

i.e. when the energy range S* becomes comparable to 
the length of quasi-energy spectrum. Numerical factor 
0.5 in ([15]) is of phenomenological origin. Indeed for data 
of Fig. 7, for L = 12, 16, 20, 24, the maxima of number 
variance lie at 9, 55, 390, 2800, whereas theoretical values 
of S* @ are 9.1, 56.9, 387.5, 2804, respectively. 
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FIG. 7. Number variance Y?(S) in log-log scale for exactly 
the same parameters (with an extra data in (b) for L = 12) 
as in previous Figure 6. 

B. Matrix elements and integrated conductance 

Knowing a complete set of eigenstates in momentum 
representation, (k\n), it is easy to compute the matrix 
elements of the current observable 
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(n\J\m) = J^(n\k)(k\m) 



(46) 



Again, one should make use of translational symmetry 
(P3|), since [S, J] = 0, to point out that matrix elements 
are non-vanishing only within a fixed quasi-momentum 
block, K n ^ K m (n\J\m) = 0. In order to obtain the 
numerical results presented in this subsection we have 
averaged over the entire Fock space (all K), except again 
for L = 24, p = 1/4 where we have averaged only over a 
block with quasi-momentum K = 1. Real part of high- 
temperature electric conductivity (g) can be (for fixed 
size L) rewritten as 

= W E \W\m)\%(^(w ~ cp m + K)) (47) 



In order to avoid ackward smoothing procedure and to 
simplify the notation we introduce scaled integrated con- 
ductivity I l {uj) 



/» = - jf a'(u)du 



another advantage of the study of direct time-evolution 
of Sec. Ill over the more common frequency-domain ap- 
proach presented here.) 




(48) 



FIG. 8. We show the integrated conductance I L (iti) for a 
cyclic chain of size L = 24 and parameter V — 2 and for 
several different values of parameter t (see legend). In the 
inset we show the same plot in semi-log scale in order to 
illustrate the zero freqency jump — the charge stiffness Dj . 



where \rj\' := min{|?/|,27r — \i]\} and 0(x) is a Heaviside 
step function. Integrated conductivity I l {(jo) is a mono- 
tonically increasing function on the frequency interval 
u> E [0, 7r], starting from the charge stiffness 



I L (0) = Dj 
and satisfying the sum-rule on the other end 



(49) 



I £ (tt) = i(J 2 ) L . (50) 



Note that the current variance can be computed |l( 



2\L 



N(L - N) 
2L{L- 1) 



lp(l_p) + 0(i) 



(51) 



In Fig. 8 we plot the integrated conductivity I L {u>) for 
different values of the kick parameter t (fixed V = 2) 
for constant size L = 24 and density p — 1/4, show- 

to non-ergodic 



ing the transition from ergodic Dj 



Dj > dynamics, consistent with results of direct time 
evolution of Sec. III. In Fig. 9 we analyze the dependence 
of I L (u) on size L for fixed values of parameters in non- 
ergodic regime t = 1, V = 2, p = 1/4. Note that for 
small frequencies cu, is rougly constant over the 

frequency interval < uo < uj^, whose width is deter- 
mined by the Thouless time of a finite system p,{L) ~ L, 
namely uo^ — 2n/L. Note that the expression for stiffness 
( fl9| ) is not completely consistent with the correct defini- 
tion (||), since the time-limit t — * oo is implicit in (^9| ) 
before TL of increasing L can be considered, whereas the 
correct order of limits is just the opposite. (This proves 




FIG. 9. We show the integrated conductance for 
a quarter filled p — 1/4 cyclic chain at t = 1, V = 2 (in 
non-ergodic regime) and for different sizes L = 24, 20, 16, 12. 
In the inset we show the same plot in semi-log scale in order 
to illustrate the convergence of the charge stiffness. Observe 
that the size of the horizontal plateau at small frequencies 
shrinks as ~ 2n/L. 

An interesting conjecture has been put forward in Refs. 
ifiof (and critically debated in 33 34|), namely that half- 
filled p = 1/2 integrable t-V model should exhibit proper- 
ties of an ideal insulator at all temperatures when V > t 
(in our notation) . Insulating behavior is characterized by 
Df = I°°(0) = and (2//?)ct'(0) = {d/dcj)I°° (0) = 0, so 
the time correlation function Cj (r) should be an oscilla- 
tory function in order to be integrated to zero. 
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We have found numerical evidence of (at least approx- 
imately) insulating behavior even in the non-integrable 
half-filled KtV model, when V > t. In Fig. 10 we demon- 
strate a double transition from (approximately) insulat- 
ing regime (example for t = 0.4, V — 1) to ideally con- 
ducting regime (example for t = V = 1) to normally 
conducting regime (example for t = 4.5, V = 1) for a 
half-filled KtV model on L = 16 sites. 




FIG. 10. Integrated conductance I L (u) for a half-filled 
(p = 1/2) cyclic chain of size L — 16 is shown for three 
different values of parameter t (and fixed parameter V = 1) 
demonstrating a double transition from insulator (here for 
t — 0.4) to ideal conductor (here for t = 1) to a normal con- 
ductor (here for t = 4.5) as the kick parameter t is increased. 
In the inset we show the same three curves in log-log scale. 



V. THIRD METHOD: ALGEBRAIC 
CONSTRUCTION OF TIME-AVERAGED 
OBSERVABLES FOR INFINITE SYSTEM 

A large amount of numerical evidence has been pre- 
sented in previous two sections in support of the Con- 
jecture put forward in Sec. I. However, all this evidence 
is based on computations on many-body systems of fi- 
nite size L, and TL has been speculated by extrapola- 
tion to 1/L = 0. One may still have doubts, that in a 
non-integrable system that is close to an integrable one 
quantum ergodicity may squeeze in for large sizes L, be- 
yond the scope of numerical observation. Therefore, as 
a complementary alternative, one would like to have a 



method of computation of time-correlators, like Da (40) 



which would directly operate with infinite systems on in- 
finite lattices L = oo. In this section we elaborate such a 
method of computation of operator valued time-average 
of an observable A in Heisenberg representation 



A = 



lim 

M-»tx 



1 



M 



2M + 1 ^ 



A(m). 



(52) 



m=—M 



The method is specially designed for kicked systems 
whoose propagators can be decomposed in several non- 



commuting parts [p5| , and will be implemented to com- 
pute time-averaged observables in infinite KtV model, in 
particular J and T", and the corresponding correlators, 
such as the charge stiffness Dj. 



A. Mathematical structures 

The first essential mathematical structure used in this 
section is the Hilbert space of pseudo-local quantum ob- 
servables. Even in the general setting we assume that 
the evolution propagator preserves the number of parti- 
cles, [U, p] = 0. So we again fix the density of particles 
p and consider observables A over a Fock (sub)space of 
quantum states with a given density parameter p. Such 
observables preserve the number of particles, [A, p] = 0. 

Let us define the scalar product of two extensive ob- 
servables A and B as 

(A\B) = lim h^B) L B = a^B) p (53) 



L — too L 



We note that ( p3[ ) has all the necessary properties of a 
scalar (inner) product: it is linear in the right factor, 
positive, and (^4|-B) = (B\A)* . Note also that averaging 
over half-filled states is in TL equivalent to the 'grand- 
canonical' average, (.) _i = (.). 

p ~2 

The observable A is called weakly local or pseudo-local, 
if ||A|| 2 := (^4 1 A) < oo. Pseudo-local observables A con- 
stitute a Hilbert space denoted by il. There is a linear 
subspace it' C it of pseudo-local observables A, such that 
[A, B] = AB — BA is pseudo-local for any pseudo-local 
B e 11. For any such A G it', the scalar product (|5^ ) 
is an invariant bilinear form with respect to the adjoint 
map (a,d A)B — [A,B], namely 

((&dA^)B\C) = (B\(&dA)C). (54) 

The second essential mathematical structure is the uni- 
tary Heisenberg-Floquet map L7 a d : il — > it, which propa- 
gates quantum observables in Heisenberg representation 
for one period of time, starting at some time rj £ [0, 1), 

U a d A(r) + m) = A{-q + m + 1) = [/ f A{rj + m)U, (55) 
{U^A\U &A B) = {A\B). (56) 

For example, for the KtV model (in spin representation 
which is, for algebraic convenience, used in this section) 
we have (19) 



U\ _i = exp(-iiriJi)exp(-iVX-ffo + ^S z )) exp(-i^tHx) 



with the potential and kinetic generator 



Hn = 



(57) 



(58) 



3+1 
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Unlike in Sec. Ill, we have here taken the time steps in the 
middle between the kicks, rj = 5, in order to fully exploit 
the time-reversal symmetry of the problem. Note that 
the time evolution of observables which are diagonal in 
momentum representation, like J and T", is not affected 
by the shift n of the origin of the stroboscopic map. The 
Floquct-Heisenberg map can be written explicitly using 
an exponential of the adjoint map 

U ad = cxp(i^tadH )exp(iVadHi)exp(i^tadH ). (59) 



with constraints 
2 



||(1 - U ad )Xr = (X\(l - U^)(l - U ad )\X) = 0, (65) 

(X\X) = const < 00, (66) 



Namely, eqs. (64-pq) imply X = a A, where a = || A\\ -1 
if X is normalized (X\X) = 1. Since constraint (f35|) is 
homogeneous the corresponding Lagrange multiplier is 
diverging. Hence, we suggest to write the constrained 
variational problem (p3J65L pq) in the compact form 



Since the density p (or magnetization A4 = p - 
spin-i formulation) is fixed, the total spin S z in eq. (|5 
generates an irrelevant phase which does not influence 
the evolution of observables. 

Time-average of observable (self-adjoint operator) A 
( [52"| ) is a solution of the fixed-point equation for the 
Floquct-Heisenberg map 



U ad A = A. 



(60) 



Time-averaging in operator space can also be written in 
terms of an orthogonal projector Pjj onto the null space 
of 1 — J7 a( j, namely 



A = P V A, P v 



The property 



lim 

Af^oo 2M 



Pu = P6 



1 M - 



m=-M 



(61) 



(62) 



is easily proved by writing a time-average limit ( pl| ) in 
an equivalent, Gaussian way 



lim j^s^X) 



X = aA, a G C 

s e (X) = ±\(X\A)\ 2 ~ lAHe-^l - e^U ad )X\\ 2 (67) 



where A is another Lagrange multiplier associated with 
the second constraint (pq). Indeed, for small e one may 
write 

s t {X) = \\(X\A)\^\\(A\A) 

_ lA( e - 2 - e - 1 )||(l-L> ad )X|| 2 + 0(e), 



so homogeneous constraint (|65|) follows automatically as 
e — > in order to make the action s t (X) regular (and 
maximal) at e = 0. Let us now show that the above 
variational problem has the correct solution (|52]). We 
differentiate the action J67| ) 



5 s e (X £ ) = (A\X e )A-±(l 



sx 



U^)(l-e-'U aA )X e 







and write a = This equation can be solved ex- 

plicitly for X e 



Pu= lim -jL- jr exp(-|(m/M) 2 )^. 



Without loss of generality we will in the following as- 
sume that observable A is traceless (A) = 0, such as J 
and T' . Note that the generalized stiffness ( f40| ) can be 
written simply as 



D A = (A\A). 



(63) 



B. Time-average from Variational Principle in 
operator space 

Scaled (or normalized) time-average of a self-adjoint 
operator A — A\ X — A/\\A\\, can be obtained from 
a variational principle in operator space, namely as an 
extremum (maximum) of an action s(X) 



SX 



s(X) = 0, 

s(X) = \{X\A)(A\X) = \\(X\At 



(64) 



X t = —(1 - e-^ad)- 1 ^ - e-U^A 
2 00 00 



ae 



E A (p) E * 

o=-oo q=\p/2\ 
,^2 



A(l-e-) 



E A (p> 



In the limit e — > 0, the last expression 
to the time-average 



X 



4a - 

lim X t = —A 

e^O A 



(68) 
is proportional 

(69) 



Evaluation 



a = (A\X) = a 



ADa 
A 



fixes the value of the Lagrange multiplier 
A = AD A . 



(70) 
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Unitarity (|56|) and invariance (|E 
important consequence 



have the following very 



(A\A) = (A\A) = (A\A). 



(71) 



Assuming that X is non-vanishing so that (X\X) > 0, 
and that X and A are proportional, one can write 



(72) 



Taking the scalar product of the last equation with A, 
one obtains a very useful expression for the stiffness 



D A 



U\x)\ 2 

(X\X) 



(73) 



C. Numerical application 

However, the maximization of the functionals (|64|, |07| ) 
in the huge infinitely dimensional operator space il is not 
convenient for practical calculation. Instead, we suggest 
to estimate the time-averaged observable A by solving 
the variational problem (|67]) in a finite-dimensional sub- 
space dJl(A) C il. (Galerkin-like approach in operator 
space). In fact, we devise a special sequence of trun- 
cated 'minimal invariant' operator spaces ...9Jl p (A) C 
9Jl p +i(A) ... C il, which in the limit p — > oo (after clo- 
sure) contain the time average A. Thus, the solutions 
X p of the variational problems (|7|) on spaces Wl p should 
converge to the proper scaled time-average X of observ- 
able A. 

Let s = {aHo + j3H\; a, (3 £ C} be 2-dim linear vector 
space spanned by the two generators of motion, H a , a = 
0,1 (|5^). Let us define the minimal invariant operator 
space containing A, as the closure of linear combinations 
of all products of adjoint generators &dH a , a — 0, 1 on 
A 



Tl(A) = \J (ads)"A 



(74) 



n=0 



Wl(A) is indeed the minimal (though infinitely- 
dimensional in general) operator space containing A with 
the invariance property 



(ad H a )m(A) =9Jt(A), a = 0,1, 

u ad m(A) = m(A), vt,v. 



(75) 



It is obvious that A £ Wl(A). We now construct the 
countable basis of the space 9Jl(A) ordered by the order 
of locality as follows: We assign an observable H q , c to an 
ordered pair of integers (q, c), order q, and code c, < 
c < 2 9 ~ 1 with q — 1 binary digits c„, c = J2n=i c n2 n_1 , 
namely 



Since not all observables H q c upto a given maximal or- 
der p,q < p, are linearly independent we perform Gram- 
Schmit orthogonalization w.r.t. the scalar product (^) 



Q c _ J Gq,c/ V {Gq,c\Gq,c)', G q ^ c ^ 0, 



0; 



G q ,c — 0, 



G„ 



(r,b)<(q,c) 
Hq.c— E G r ,b(G ri b\H q 



(77) 



The nonzero observables G g>c form the orthonormal basis 
of VJl(A). Note that observables G 9iC are strictly local 
operators of order q: in case of spin representation of 
KtV model they are represented as expansions 



G 



q,c 



E az 



(78) 



in terms of spatially homogeneous finite products of field 
operators 



Z x 



= E 



where Sk G {0, +, — , z} and cr^ = 1. The (average) num- 
ber of nonzero terms in expansions (f78j) was found to 
grow exponentially at approximately the same rate for 
both observables under study, J and T', namely as 



WqT 



± 0} w 0.5 x 2.55 9 



which may be further reduced by a factor 2, or even by a 
factor 4 if p = ^ , using the symmetries V (J2J) , and 71 ( |25| ) 
(the latter may be used if p = |). Note that the entire 
linear space VJl(A) has the same symmetry properties as 
observable A, for example the space Wl(J) and 97t(T') 
belong to a negative and positive parity symmetry class, 
respectively, w.r.t. parity operation (p4|), 

Let us now define a sequence of truncated minimal in- 
variant operator spaces containing A, 



p-i 



m p (A) = |J(ads)"A, p=l,2. 



(79) 



71=0 



with dimensions d p (A) := dimWl p (A). Linear space 
dJl p (A) contains operators derived from A by composi- 
tion of generators ad H a up to order p. Due to binary 
code construction (|7^) we have a strict upper bound on 
the growth of dimensions of the truncated spaces 9Jl p (A), 



d p (A) < 2 p - 



(80) 



however actual dimensions may grow considerably slower 
(due to systematic linear dependences among -ffp,f>), 
namely for A — J and A — T' we find by computer 
algebra up to p = 14th order (see Tab.l) 



H q , c = (adiJ C9 _ 1 )(adi? C5 _ 2 ) • • • (adff Cl )A 



(76) 



dp (J) « 1.825 p -\ d p (T') « 1.68 p " 



IV. 
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Let Hp. ai a = 0,1, denote real and symmetric (Her- 
mitean in general) matrices of linar maps &dH a on 
Tlp(A) with images orthogonally projected back to 
Tlp(A). It follows from the construction that they have 
(generally) a block-banded structure where the blocks 
correspond to observables with a fixed order q: namely 

{G q , c \adH a \G ql , c ,) = 0, if \q-q'\^l. (82) 

The truncated adjoint maps, H p a , have nontrivial null 
spaces 

%,, a {A) = {B e Ttp(A); [H a , B] G Wl p+1 (A) - Wl p (A)}, 

(83) 

with dimensions d Pt a(A) := dim91 p , Q (A) which increase 
approximately with the same exponent as d p (A) ( |8l"| ) (see 
Tab.l). 

By means of truncated adjoint maps H p . a we construct 
an approximate Floquet-Heiseberg matrix U p , which is 
a dp (A)-dim unitary matrix over the truncated space 

m p {A) 

Up = exp(ii<H Pi i) exp(zFHp i0 ) exp(iiiH Pi i) (84) 

Now we are ready to solve the variational problem ( |64j - 
|67]) in the truncated space 9Jl p (A). We note an impor- 
tant 'experimental' observation (whose theoretical under- 
standing is stil lacking) namely that the matrix 1 — U p 
possesses a high-dimensional null-space 

Vlp (A) = {Be m p (A);V p B = B}, 

whose dimension d p (A) :— dim^tp 1 (A) is, for odd p, in- 
dependent of parameters t, V and equal to the dimension 
of the null space of H Pi i, 

d%_ 1 (A) = <h l - ltl (A). (85) 

Note also that for odd order of truncation p = 21 — 1, 
the elements of null space B E 91^ (A) are spanned 
by combinations of even powers of generators only, i.e. 
(B\G2i, c ) = 0, which is due to time-symmetric construc- 
tion (rj = i) of the evolution operator C/ ad (^9|). 

The scalar products (|5^) for different values of the den- 
sity p are non-degenerate with respect to each other, and 
therefore the dimensions of various linear (sub)spaces, 
d p (A), d Pta (A), dp (A), (see Tab.l) do not depend on the 
density parameter p. 

The constraint ( |65| ) is now equivalent to restricting the 
variation j6~i| ) to the sub-space 9t p (A). Hence, the 'trun- 
cated' scaled time-averaged observable X p is a maximum 
of the quadratic form (X p \A)(A\X p ) on Vl u (A) with a 
normalization constraint (X p \X p ) = 1. In other words, 
if F n ,n — 1, . . . , d :— d p (A) is an orthonormal basis of 
the null-space yt p (A) and if (xi, . . . , Xd) is a normalized 
eigenvector of the (positive definite) d x d matrix eigen- 
value problem 

J2{Fm\A){A\F n )x n = fx m 

n 



with the maximal eigenvalue /, then 

Xp=Y, F nXn (86) 
a 

is a solution of the variational problem (|64|-p7|) in the 
truncated space 9Jl p (A). In the limit p — > oo we expect 
to recover an exact time-average 

lim X p =X = || A (87) 

p— >oo 

However, if the system is ergodic, time average should be 
zero A = (note that (A) — 0), so the (normalizable) 
limit ( |S7| ) should not exist. In order to inspect the con- 
vergence of X p in the Hilbcrt space of observables 11, we 
define a relative norm N q (X) w.r.t. order q 

N q (X) = Y,\(X\G q , c )\ 2 . (88) 

C 

Since 

oo 

||X|| 2 = 5> g (X), (89) 

9=0 

the inspection of the convergence of the sum on RHS of 
( |S9| ) would give us indication of convergence of X p flslf ) 
and thus non-ergodicity of the problem. As a second 
criterion of convergence of X we study stability of X p , or 
of the relative norms N q (X p ),q < p, w.r.t. variation of 
the truncation order p. 



p 


d p (J) 




d p ,i(J) 


d P (T') 


dpfi(T') 


dp AT 1 ) 


1 


1 


1 


1 


1 
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2 
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2 
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12 
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10 
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21 
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15 
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38 


16 


12 


25 


9 


10 


8 


69 


27 


21 


40 


12 


7 


9 


126 


48 


38 


66 


22 


21 


10 


230 


84 


68 


107 


33 


22 


11 


419 


153 


123 


178 


60 


51 


12 


763 


273 


223 


293 


91 


66 


13 


1393 


493 


409 


494 


162 


137 


14 








831 


257 


202 



TABLE I. Dimensions of the truncated minimal invariant 
spaces and of the null-spaces of truncated adjoint maps for 
different orders of truncation p. 
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In Figs. 11,12 we show the relative norms N q (X p ) of 
the normalized time-average of the current J (Fig. 11) 
and kinetic energy T' (Fig. 12) for several different orders 
of truncation p (up to p — 11 for J and up to p = 13 for 
T'). We note that for both observables, X p is quite stable 
against variation of p, for t < 1.4, and in the same param- 
eter range, the coefficients N q (X p ) seem to be summable. 
The stiffness Da ( |73] ) may be re-written in terms of rel- 
ative norms as 



D A = (A\A) 



Ni(X) 



(90) 



In this regime where X is convergent in it, N q {X p ) are 
good approximation of N q {X) for q < p (appart from 
a constant renormalization prefactor which cancels out 
from expressions like (|90|)), and we may write a good 
estimate for the upper bound on the stiffness 



U A — 



\{A\X V 
> (A\A) 



2 a (A\A) 
N X (X) 



ET =1 N q (X) 



Y,UN q {X) 
Df. 



(91) 



However, we would like 
tions of the stiffness D c , 



to have accurate approxima- 
y A itself rather than the accu- 
rate upper bounds, so we extrapolate the relative norms 
N q (X) w Ng(Xp) to orders higher than the order of trun- 
cation, q > p, in expression for the stiffness ( |90| ) by fitting 
the tail of N q=2 i+i (X p ) at three points: q = p — 4, p — 2, p 
(note that N q= 2i(X p ) — 0) with exponential ansatz 
N q= 2i+i{X p ) cx exp(— sq). Since the actual rate of con- 
vergence of N q (X) — > 0, as q — > oo is probably slower 
than exponential (see Figs. 11, 12), the stiffness extrapo- 
lated in this way (eq. (|90|)), D A , is probably still slightly 
overestimated. In Fig. 13 we show the dependence of the 
extrapolated charge stiffness D e j on the parameter t (and 
for fixed V = 2) through the critical range t ~ t c ra 1.45, 
and compare it with the charge stiffness as computed 
from direct diagonalization of the finite KtV chains of 
sizes L — 24, 20, 16. When approaching the critical point 
t c , the fitted slope s linearly decreases to zero. For 
larger values of parameters, t > t c (V), X p is not sta- 
ble against variation of p and the partial sums of relative 
norms N q (X) are not converging. Therefore, A = and 
— 0, and the system is quantum ergodic. In Fig. 14 
we plot a full (t, V) phase diagram of the (extrapolated) 
charge stiffness Dj. It is clear that this last method, since 
in operates with an infinite system, gives the most reli- 
able results on the critical regions of transition between 
dynamical phases. However, no other dynamical infor- 
mation on correlation functions is obtained other than 
their time-averages, so within the present method we can- 
not make any claims on the stronger ergodic property of 
quantum mixing. 

Although dynamical behavior of observables may de- 
pend on the symmetry class of observables w.r.t., say, 
parity operation (pi|), we have found very similar ergodic, 



non-ergodic, and critical regions, for the two examples of 
opposite parity observables, J and T", that have been 
studied. However, we should note that dynamical be- 
havior may also depend on the other (continuous) con- 
served quantities, such as the density p. Our results for 
other values of density p indicate that the transition re- 
gion between ergodic and non-ergodic dynamics moves 
to slightly smaller values of parameter t as p approaches 
1/2. Due to paricle-hole transformation ( p5[ ) the dynam- 
ics for p = p' is equivalent to dynamics for p = 1 — p'. 




FIG. 11. The logarithm of relative norm N q (X) of the nor- 
malized time-averaged current X — J/\\J\\ in a quarter-filled 
(p = 1/4) infinite KtV model is plotted against (odd) order 
q = 21 — 1, for a square mesh of parameters t and V . The 
three different curves on each graph, thick, medium, and thin, 
refer to three different orders of truncation of operator spaces, 
p = 11, p = 9, and p = 7, respectively. 

We should note that in a recent paper |36| a very sim- 
ilar algebraic approach has been used in order to com- 
pute the pseudo-local quantum invariants of motion. In 
the regime of non-ergodic dynamics one or two converged 
pseudo-local invariants of motion were found, whereas in 
the regime of ergodic dynamics, no non-trivial invariants 
of motion were found. Then by using a formula (|9|) of 
Mazur jl3| and Suzuki the time averaged correla- 
tion of kinetic energy Dt has been computed by means of 
an expansion in terms of pseudo-local invariants, giving 
the results which are in good agreement with the results 
of direct calculations on finite systems. We believe that 
our variational approach in the space of observables pre- 
sented here is (in general, possibly non-integrable case) 
an improvement of the Mazur-Suzuki approach jl3 14 
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to the calculation of time- averaged correlators. Within 
Mazur formula (|J) one is typically able to write only in- 
equality (lower bound on stiffness) since the set of known 
local invariants of motion may be incomplete. 




1 3 5 1 9 13 q 

V-0.4 V-1.0 V-1.6 V-2.2 V-2.8 



FIG. 12. Same as in previous Fig. 11 for the normalized 
time-averaged kinetic energy X — T'/\\T'\\, except for larger 
truncation orders, p = 13 (thick curves), p = 11 (medium 
curves), p — 7 (thin curves). 




FIG. 13. Charge stiffness Dj vs. kick parameter t and 
constant parameter V — 2 for quarter-filled chain p = 1/4. 
Different curves refer to different system sizes L = 24, 20, 16, 
while points represent the infinite-size stiffness (j?^) based on 
extrapolated algebraic time-averaged current invariant of mo- 
tion. The truncation order is p = 11. In the inset we show 
the logarithmic slope s of the falloff of the relative norms 
N q (J) oc exp(— sq) at q ~ p. 



Current Invariant 




FIG. 14. (t, V) phase diagram of the charge stiffness Dj 
for quarter-filled (p — 1/4) infinite KtV model, as deduced 
from extrapolated time-averaged current J. 



VI. CONCLUSIONS AND DISCUSSION 

In this paper we have presented three complementary 
(mainly numerical and computer-algebraic) approaches 
to the dynamics of non-integrable quantum many-body 
systems in thermodynamic limit (TL), demonstrated and 
studied in a kicked t-V model of spinless fermions. We 
have been primarily interested in the structural stability 
of non-ergodic quantum motion and the transition from 
non-ergodic/non-mixing to ergodic/mixing dynamics in 
TL. 

The first approach that we have used is a direct time- 
evolution of a finite quantum system (which may be in 
the present model performed very efficiently by means 
of the so-called Fermionic Fast Fourier Transformation) 
and computation of time-correlation functions of generic 
quantum observables. The size L of the system is sys- 
tematically increased and TL is speculated based on ex- 
trapolation to 1/L = 0. For sufficiently large values of 
kick parameters, we have found quantum mixing and ex- 
ponential decay of time correlation functions, while for 
smaller, intermediate values (~ 1) of kick parameters, we 
have found non-mixing quantum motion characterized by 
saturating, non- vanishing time-correlation functions. 

Our second approach is a direct diagonalization of the 
stationary quantum problem of finite size and calculation 
of dynamical properties, such as charge stiffness, conduc- 
tivity, etc., in frequency domain. Also traditional quan- 
tum signatures of chaos, such as level statistics, has been 
inspected and shown to correspond with the dynamical 
behavior. This approach is less computationally efficient 
in case of the present model than the first one. 

In the third approach, which is fully complementary 
to the other two, we propose an algebraic method for 
computation of time-averaged observables of an infinite 



l(i 



system. Thus we can make most precise statements on 
quantum-ergodicity of an infinite system, which are in 
complete agreement with the extrapolated results of cal- 
culations on finite systems. 

The above results are claimed to be the evidence for 
the validity of the Conjecture (Sec. I), namely that the 
generic quantum-many-body system in TL may not be 
quantum-ergodic (or mixing) if it is sufficiently close to 
an integrable system in parameter space. Recent nu- 
merical results on transport in extended (non-integrable) 
Hubbard model |33| are compatible with the above Con- 
jecture. The transition between non-ergodic and ergodic- 
dynamics when the external parameters are increased 
has the properties of (dynamical) phase transition and 
should be further studied theoretically. First such at- 
tempt has been undertaken in j57|, where a discretized 
non-integrable quantum field model (in the continuum 
limit) has been mapped on a quantum chaotic model of 
a single particle on a 2-dim torus (in the quasi-classical 
limit), and the transition from non-ergodic/non-mixing 
dynamics to ergodic/mixing dynamics of the quantum 
field model has been identified with the stochastic tran- 
sition from regular to chaotic motion. 

We have also given a clear evdience on the non-trivial 
existence of mixing quantum motion in KtV model in TL 
with exponentially decreasing time-correlation functions, 
provided the external (kick) parameters are large enough 
(above the critical values). Such quantum mixing behav- 
ior may be a source of truly chaotic and macroscopically 
irreversible quantum motion of many-body systems p8[ . 
Note that macroscopic irreversibility as a consequence 
of non-dissipative but strongly non-integrable quantum 
many-bod y d ynamics has been recently observed exper- 
imentally |39| . 

One might doubtfully argue that our quite surpris- 
ing finding on structurally stable non-ergodic quantum 
motion in TL (formulated as Conjecture in Sec. I) may 
be just another peculiarity of Physics in One-dimension, 
and as such should not be expected to hold in inter- 
acting quantum systems in more than one spatial di- 
mension. Being aware of this fear we have straightfor- 
wardly extended our KtV model ( |To| , p~T| ) to a rectangu- 
lar periodic Li x I2 lattice in two spatial dimensions, 
with isotropic hopping in two orthogonal directions and 
5— kicked isotropic nearest neighbour interaction. An ef- 
ficient direct time evolution of the 2-dim KtV model has 
been implemented analogously along the lines described 
in Sec. Ill, and its time correlation functions have been 
computed accordingly, however due to greater computa- 
tional complexity only for relatively small lattices of sizes 
up to 6 x 5. We should stress that we were again able to 
identify quite clearly the two regimes of quantum motion 
which have been roughly stable against the variation of 
the lattice size, namely: (i) a quantum mixing regime for 
sufficiently large t and V, and more important, (ii) quan- 
tum non-ergodic and non- mixing regime for \t\ ~ |V| ~ | 
(or smaller), although the system is not known to be an- 
alytically integrable in the limit £ — ^ 0, — ^ 0. This 



result (whose details will be published elsewhere) is a 
small piece of numerical evidence for the validity of the 
Conjecture in two spatial dimensions. 

Discussions with Prof. P. Prelovsek, and the financial 
support by the Ministry of Science and Technology of R 
Slovenia are gratefully acknowledged. 
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